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Advanced Small Perturbation Potential Flow Theory for 
Unsteady Aerodynamic and Aeroelastic Analyses 

John T. Batina' |, 

NASA Langley Research Center 
Hampton, Virginia 23681 

Abstract 

An advanced small perturbation (ASP) potential flow theory has been developed to 
improve upon the classical transonic small perturbation (TSP) theories that have been 
used in various computer codes. These computer codes are typically used for unsteady 
aerodynamic and aeroelastic analyses in the nonlinear transonic flight regime. The codes 
exploit the simplicity of stationary Cartesian meshes with the movement or deformation 
of the configuration under consideration incorporated into the solution algorithm through 
a planar surface boundary condition. The new ASP theory was developed methodically 
by first determining the essential elements required to produce full-potential-like 
solutions with a small perturbation approach on the requisite Cartesian grid. This level of 
accuracy required a higher-order streamwise mass flux and a mass conserving surface 
boundary condition. The ASP theory was further developed by determining the essential 
elements required to produce results that agreed well with Euler solutions. This level of 
accuracy required mass conserving entropy and vorticity effects, and second-order terms 
in the trailing wake boundary condition. Finally, an integral boundary layer procedure, 
applicable to both attached and shock-induced separated flows, was incorporated for 
viscous effects. The resulting ASP potential flow theory, including entropy, vorticity, 
and viscous effects, is shown to be mathematically more appropriate and computationally 
more accurate than the classical TSP theories. The formulaic details of the ASP theory 
are described fully and the improvements are demonstrated through careful comparisons 
with accepted alternative results and experimental data. The new theory has been used as 
the basis for a new computer code called ASP3D (Advanced Small Perturbation - 3D), 
which also is briefly described with representative results. 

Introduction 

Classical transonic small perturbation (TSP) potential flow theories have been used 
over the past three decades or so in the development of computer programs for unsteady 
aerodynamic and aeroelastic analyses. 1 ' 3 Notable efforts in this regard are the 
XTRAN3S 4 code developed by Boeing under United States Air Force sponsorship and 
the CAP-TSD 5 code developed at the NASA Langley Research Center. These codes 
were developed primarily for the investigation of aeroelastic phenomena in the nonlinear 
transonic flight regime such as flutter and divergence. The codes exploit the simplicity of 
stationary Cartesian meshes with the movement or deformation of the configuration 
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under consideration incorporated into the solution algorithm through the surface 
boundary condition. This approach is in contrast to the more sophisticated modeling 
incorporated into advanced computer codes developed to solve the full potential, Euler, 
and Navier-Stokes equations, wherein the body-fitted mesh deforms locally to the 
instantaneous shape or position of the configuration. 6 However, the latter approach is 
computationally more expensive, and when added to the generally much higher cost of 
solving the higher-order equations (Euler and Navier-Stokes) themselves, it is cost 
prohibitive for routine aeroelastic application. Therefore, the small perturbation approach 
is preferred because of the multitude of cases that typically need to be considered for 
aeroelasticity. 

Over the years research has been conducted by several organizations including the 
NASA Langley Research Center to assess the general applicability and accuracy of the 
TSP codes in numerous unsteady aerodynamic and aeroelastic applications. 7 ' 10 These 
applications have revealed various inaccuracies inherent in the underlying small 
perturbation theories 11 including accuracy limitations near the leading edge and wing tip 
regions, 7 shock waves that are inaccurately captured in terms of strength or location, 10 
and convergence difficulties attributable to certain computational and mathematical 
complications such as non-uniqueness of the potential flow equation. 12 Consequently, the 
author conducted research over the last decade to examine the applicability of the small 
perturbation concept in general and the accuracy of the CAP-TSD code in specific. The 
author is the principal developer 13 ' 15 of the CAP-TSD code and is therefore in a unique 
position to make a critical assessment of the methodology contained therein. The 
objective of the present effort was to identify the sources of the inaccuracies and 
determine if improvements could be made to alleviate or eliminate them. 

The first step in this effort was to determine the accuracy of the CAP-TSD code for 
the in viscid flow about a basic two-dimensional airfoil configuration. The NACA 0012 
airfoil was selected because the shape is analytically defined and calculations could be 
performed first for non-lifting cases to eliminate any potential inaccuracies due to the 
modeling of the trailing wake. There are standard cases 16 for the NACA 0012 airfoil 
whereby accurate solutions, obtained using computer codes that solve the full potential 
(FP), Euler, and Navier-Stokes equations, have been published. These solutions were 
used for comparison purposes. The inviscid CAP-TSD calculations differed substantially 
from the accepted FP and Euler solutions. For transonic cases, the shockwaves were 
poorly predicted in terms of strength and location. This was the case even when the 
shock-induced entropy and vorticity effects 15 were included in the CAP-TSD calculation 
and comparisons were made with accepted Euler solutions. For subcritical cases where 
such effects are unnecessary, the pressure levels were poorly predicted in comparison 
with accepted full-potential solutions even when the complete pressure coefficient 
formula was selected within the code. This suggests that the underlying TSP potential 
flow theory in CAP-TSD is deficient in some respect, and the inclusion of entropy, 
vorticity, and even viscous effects is then only an academic exercise. 

The next step in this effort was to determine the source of the inaccuracies and 
develop modifications to hopefully improve the accuracy of the resulting capability. As a 
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result, a new so-called advanced small perturbation (ASP) theory was developed 
methodically by first determining the essential elements required to produce FP-like 
solutions with a small perturbation approach on the requisite Cartesian grid. This level of 
accuracy was obtained by using a higher-order streamwise mass flux and a mass 
conserving surface boundary condition. Subsequent applications with these 
improvements showed very good agreement with all of the full potential cases 
considered. This was true even for cases at freestream Mach numbers as low as 0.3 and 
angles of attack as high as ten degrees, conditions that are normally considered to be 
outside the range of applicability of classical small perturbation theories. It was then 
discovered that when the CAP-TSD entropy/vorticity modeling 15 was added to the ASP 
modifications, the resulting shock waves did not have the correct strength and location in 
comparison with accepted Euler solutions. Further research indicated that an essential 
element to ASP theory (and the correct approach for any small perturbation theory for 
that matter) to produce Euler-like results is a mass conserving implementation of the 
entropy effects. In CAP-TSD the Rankine-Hugoniot (R-H) shock jump relation 17 is used 
to calculate the entropy jump across the shock. However, this approach does not 
conserve mass because it is not consistent with the small perturbation governing 
equation. (The R-H calculation is mass conserving only if the governing equation is the 
FP equation.) Thus a mass conserving calculation of the entropy jump across of the 
shock wave was derived based on the higher-order streamwise mass flux of the ASP 
theory. Second-order terms in the trailing wake boundary condition were also found to 
be not insignificant for lifting cases and thus were incorporated. Subsequent calculations 
with these modifications showed very good agreement with all of the Euler cases 
considered. 

The final step of the effort was to implement viscous effects by coupling an integral 
boundary layer procedure to the outer inviscid flow calculation. The dissipation integral 
method 18 ' 20 was implemented to model attached or shock-induced separated boundary 
layers. The capability involves solving the unsteady boundary layer equations 
simultaneously with the outer potential flow solution, so that no interaction law coupling 
the inner and outer solutions is required. The combined solution procedure involves an 
implicit block tridiagonal inversion for all of the cells along the surface and its trailing 
wake. Unlike the CAP-TSD viscous capability, 21 exact formulas for edge quantities and 
exact boundary conditions are used along surfaces and wakes. Smoothing of edge 
quantities and limiters are not required for stability, and no arbitrary or free parameters 
are necessary to tune the procedure. Calculations performed with the complete ASP 
capability including entropy, vorticity, and viscous effects showed good agreement with 
experimental data. 

Formulaic details of the newly developed ASP potential flow theory including 
entropy, vorticity, and viscous effects are described fully. The ASP concepts are 
demonstrated through careful comparisons with accepted alternative results and 
experimental data for the NACA 64A410, NACA 0012, and RAE 2822 airfoils. Most of 
the results are for steady flow applications to prove the ASP concept. A separate 
publication will report unsteady calculations. The ASP theory is shown to be 
mathematically more appropriate and computationally more accurate than the classical 
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TSP theories. Consequently, the new theory has been used as the basis for a new 
computer code called ASP3D (Advanced Small Perturbation - 3D), which involves either 
AF1- or AF2-type approximate factorization algorithms and a FAS (full approximation 
scheme) multigrid procedure for the solution of the ASP potential flow equation. The 
ASP3D code can treat aircraft configurations involving multiple lifting surfaces with 
leading and trailing edge control surfaces and a fuselage. The new code is briefly 
described and some representative results are presented. Detailed descriptions of the 
ASP3D capability and comprehensive results and comparisons are reported in a 
companion publication. 22 

Advanced Small Perturbation Theory 

In this section the essential elements of the ASP theory are described, along with 
mathematical and computational comparisons with various TSP theories, selected 
alternative results, and experimental data. 

Governing Equation 

It is first instructive to examine in detail the general small perturbation equation that 
is presented here in two dimensions for simplicity as 

* + 4 . + *_0 

dt dx dz 

The fluxes are defined as 

/ 0 =-A^-Z% 
fi = C + D(f) x + E(f>l + F(j>l 

f3 = t 


The constant coefficients are defined as 

A = Ml B = 2 Ml C = 1 D = 1 - Ml 

The remaining coefficients E and F are somewhat arbitrary depending upon the 
assumptions used in deriving the governing equation. When the equation is derived to 
match a small perturbation form of the shock condition the so-called NASA Ames 
coefficients result, given by 

£ = _I (y+ l) M 2 F=0 

The so-called NLR coefficients , derived by taking a small perturbation approximation of 
a series expansion of the density, are defined by 
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F = 0 


E = ~\[ 3-(2- r )Ml]Ml 

And the newly defined ASP ( Advanced Small Perturbation) coefficients, derived by 
matching the exact sonic and stagnation conditions, are defined by 

£ = -I (y+ l) M 2 F= _I (y + l) M 2 

2 o 


An alternative streamwise flux f was derived by an asymptotic expansion of the 
Euler equations by Williams. 23 This alternative flux is defined by 


where 


V = 


/l = (/ + l)Af»[l + (<l> x ) sonic] 
<t>x 
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By design, Williams’ flux has the exact sonic velocity and the correct shock jump 
condition. Therefore it is a good alternative to the TSP mass flux evaluated using the 
NLR or NASA Ames coefficients for cases involving shock waves. However as 
discussed below, Williams’ flux is inaccurate near stagnation, which results in numerical 
difficulties near the leading edge. The formulas for the streamwise mass fluxes and 
associated coefficients and parameters from the various small perturbation theories are 
tabulated in the upper part of Table 1 for easy reference and direct comparison. 

To compare and contrast the various small perturbation formulations graphically, 
Figure 1 shows flux quantities as functions of streamwise velocity (u = 1 + (p x ) for the 
ASP theory, the TSP theories with NLR and NASA Ames coefficients, and Williams’ 
asymptotic expansion. The flux quantities were computed at a freestream Mach number 
of M v = 0.72. Figure 1(a) shows the streamwise mass flux f\ and Figure 1(b) shows the 
derivative of the streamwise mass flux df/d<p x . The derivative of the mass flux is 
important because it is related to the product of the local wave speeds given by 

= M 2 (a - u)(a + u) = M 2 (a 2 - u 2 ) 

d <Px 

where a is the speed of sound. The sign of the derivative df/d<p x indicates whether the 
local flow is subsonic (positive because \ u\ <a) or supersonic (negative because | u\ > a). 
And hence, the sonic point is the velocity by which df/d<p x = 0. 

At u = 1, which corresponds to the undisturbed or freestream flow of (p x = 0, the four 
small perturbation theories have identical values of the flux and its derivative, as 
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expected. Near u = 1, which corresponds to small perturbations | (p\ « 1 where classical 
small perturbation theory is mathematically valid, the flux quantities are similar between 
the various theories with the greatest similarity occurring between the TSP theories with 
the NLR and NASA Ames coefficients. For larger perturbations corresponding to 
approximately w<0.8orw>1.2, there are larger differences between the flux quantities 
of the small perturbation theories, especially for slower speeds u -* 0 corresponding to 
stagnation, and even more so for reverse flows u < 0. Formulas for the derivatives of the 
mass fluxes and associated sonic velocities from the various small perturbation theories 
are tabulated in the lower part of Table 1 for easy reference and direct comparison. Table 
1 emphasizes that the sonic velocity is exact for the ASP theory and Williams’ 
asymptotic expansion theory. 

Of the various small perturbation formulations shown in Figure 1, the ASP theory is 
the only formulation that has the correct form for the mass flux and its derivative across 
the entire velocity range. The correct form for the mass flux is an asymmetric function 
about u = 0 (stagnation) because the mass flux physically should have a similar 
magnitude but a negative sign for reverse flow as shown by the ASP theory curve in 
Figure 1(a). The TSP theories with the NLR or NASA Ames coefficients and Williams’ 
formulation do not satisfy the above property, and hence, they are not applicable for cases 
involving reverse flow. Furthermore, they are inaccurate near stagnation, because those 
formulations have negative values for the mass flux at small positive values of u, with the 
greatest deviation from zero given by the Williams’ flux. 

The correct form for the derivative of the mass flux is a symmetric function 
about u = 0 (stagnation), as shown by the ASP curve in Figure 1(b). This is because the 
derivative should be positive for subsonic flow, corresponding to | u\ < a, which includes 
reverse subsonic flow, and the derivative should be negative for supersonic flow, 
corresponding to | u\ > a, including reverse supersonic flow. When the derivative of the 
mass flux is positive, the governing equation is of elliptic type, correctly describing 
subsonic flow. When the derivative is negative, the governing equation is hyperbolic, 
physically corresponding to supersonic flow. For the ASP theory this is true, 
independent of the direction of the flow. The other theories not only have the wrong 
form but they are inaccurate near stagnation and inaccurate for stronger supersonic flows. 
In fact the worst formulation in this regard is the William’s flux, which explains why 
applications of CAP-TSD using the William’s flux (default formulation) encounter 
convergence difficulties, especially on finer meshes. The formulas for the speed of 
sound, velocity, and stagnation quantities, derived from the derivatives of the mass flux 
for the various small perturbation theories, are included in Table 2. The entries for 
Williams’ theory are marked “not available” because the quantities cannot be derived 
easily from the asymptotic mass flux. All of the quantities tabulated for the ASP theory 
are exact because the derivative of the mass flux dfjd(p x is exact. Hence, applications of 
the ASP theory are expected to be accurate for conditions ranging from stagnation to 
sonic to strong supersonic, including any reverse flow that might occur. 

Figure 2 shows flux quantities as functions of streamwise velocity (u = 1 + <p x ) for 
the ASP theory computed for freestream Mach numbers ranging from M x = 0.0 to 1.2. 
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Figure 2(a) shows the streamwise mass flu x/i and Figure 2(b) shows the derivative of the 
streamwise mass flux df\/d(p x . Figure 2(a) emphasizes that the ASP mass flux is an 
asymmetric function about u = 0 (stagnation) and the M x = 0.0 result is a line with a slope 
of unity. If desired the positive offset in/i can be removed by replacing the constant C = 
1 in the ASP mass flux by the constant C = D - E + F, although this modification was not 
investigated. Of course the derivative is independent of the constant C. Figure 2(b) 
emphasizes that the derivative of the ASP flux is a symmetric function about u = 0 
(stagnation) at any freestream Mach number. 

Figure 3 shows an expanded view of the derivative of the mass flux as a function of 
the streamwise velocity for the M rx = 0.72 case of Figure 1 to illustrate the differences in 
sonic points predicted by the various small perturbation formulations. The sonic point is 
the velocity by which dfjd(p x = 0. As discussed before, the ASP and Williams’ flux 
formulations possess the exact sonic point by design as shown in Figure 3. The TSP 
theories with the NLR and NASA Ames coefficients predict values of the sonic velocity 
that are too high, with the highest value of the sonic velocity corresponding to the NASA 
Ames formulation. The TSP theory with the NASA Ames coefficients will consequently 
produce transonic solutions with the smallest supersonic zones and weakest shock waves 
in comparison with the other formulations. This is due to the flow remaining subsonic 
until the local velocity exceeds the higher predicted value of the sonic velocity. The TSP 
theory with the NLR coefficients will generally produce transonic solutions with larger 
supersonic regions and stronger shocks in comparison with the TSP-NASA Ames theory, 
but the shock waves will not be as strong as those predicted correctly by either the ASP 
or Williams’ formulations. 

Flow calculations were performed with the small perturbation theories using a 257 x 
129 Cartesian finite volume mesh that had a fifty chord extent as shown in Figure 4. The 
spacing at the leading and trailing edges is Ax = 0.005c, as shown in Figures 5(a) and 
5(b), respectively, which is generally considered small for transonic small perturbation 
applications. All of the solutions were converged at least seven orders of magnitude in 
the reduction of the L 2 -norm of the residual. The effects of the streamwise flux on the 
pressure coefficient distribution for a typical transonic case are shown in Figure 6 for the 
ASP theory and the TSP theories with the NLR and NASA Ames coefficients. Results 
were not obtained using Williams’ asymptotic flux because of numerical difficulties in 
the leading edge region that inhibited convergence. The case considered is that of the 
NACA 64A410 airfoil at an angle of attack of a = 0° and M v = 0.72, the same freestream 
Mach number used in Figures 1 and 3. As shown in Figure 6, the solution obtained using 
the TSP theory with the NASA Ames coefficients indicates that there is a smaller 
supersonic region and a weaker shock wave on the upper surface of the airfoil in 
comparison with the other two solutions as predicted by the theoretical analysis discussed 
earlier. The solution obtained using the TSP theory with the NLR coefficients has a 
slightly stronger shock wave, but the strongest shock wave is predicted by the ASP 
theory, also consistent with the theoretical analysis. The ASP pressure distribution is in 
very good agreement with a conservative full-potential (FP) calculation shown in Figure 
7, reported by Jameson 24 for this case. A comparison of the ASP (Figure 6) and FP 
(Figure 7) pressure distributions demonstrates that the strength and position of the shock 
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wave near approximately 63% chord on the upper surface of the airfoil are predicted 
accurately by the ASP formulation. Hence, the ASP mass flux is an essential element in 
obtaining FP-like solutions with a small perturbation formulation. 

Inviscid Irrotational Surface Boundary Conditions 

The CAP-TSD code and other small perturbation codes typically use the lowest order 
classical small disturbance boundary condition to impose surface tangency given by 

(p z = b x -a 

where 

b(x ) - xa = 0 

represents the surface. With this boundary condition the velocity and mass flux vectors 
do not generally coincide, and neither vector is tangent to the surface. 

Use of this boundary condition generally results in shock waves that are too weak 
and located forward in comparison with full potential results. Furthermore such solutions 
can be shown to be dependent on the local mesh density. Specifically, if the mesh is 
changed, such as using a finer or coarser mesh, the shock strength and location will be 
materially changed. 

Requiring that the velocity vector be tangent to the surface results in 


<Pz =( 1+ 0x)(^-«) 


which may be referred to as a “full potential” surface boundary condition, similar to what 
is imposed in full potential codes albeit on a body-fitted mesh. However this boundary 
condition is not consistent with the governing equation and generally leads to solutions 
involving shock waves that are too strong and located too far aft. 

Requiring that the mass flux be tangent to the surface results in 

<p z = f\(b x -a) 

where f, is the streamwise mass flux. This boundary condition is consistent with the 
governing equation (mass consistent) and its use produces solutions that are relatively 
mesh independent, meaning that the local mesh density does not appreciably affect the 
shock strength and location. 

The ASP capability uses a higher-order mass conserving boundary condition to 
impose surface tangency given by 

<P Z = ~( b x ~ a ) 
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g = 1 + H(p x + and H = -(y-l)Ml 

Note thatg = Mia 2 . This boundary condition was found to be more accurate than the 
simpler one above, although the inclusion of the complete vertical flux g in the governing 
equation had a negligible effect on the solution. 

The effects of surface boundary condition (BC) on pressure coefficient distribution 
are shown in Figure 8 for three forms of the surface BC including the lowest order form 
used in CAP-TSD labeled “slopes”, the full potential BC where the slopes are scaled by 
the streamwise velocity (u = 1 + <pj, and the mass consistent BC used in the ASP 
formulation. The calculations were performed for the same NACA 64A410 case 
presented in Figure 6 and the ASP streamwise flux was used throughout. Figure 8 shows 
that the use of the lowest order “slopes” BC produces a solution that has a shock wave 
that is too weak located forward of the correct position. Although not shown here, 
solutions obtained with the lowest order BC tend to be mesh-dependent because the BC is 
not mass conserving. Use of the FP BC produces a shock wave that is too strong and 
located too far aft. In contrast, the mass consistent BC (with the ASP mass flux) leads to 
the correct shock strength and location, as compared with the Jameson 24 FP solution of 
Figure 9 (same as Figure 7; repeated for direct comparison). Hence, the mass consistent 
surface boundary condition is an essential element in obtaining FP-like solutions with a 
small perturbation formulation. 

Inviscid Irrotational Wake Boundary Condition 

The wake is modeled in all of the small perturbation theories as a planar shear layer 
emanating from the trailing edge of the wing. The circulation in the wake is first 
computed at the trailing edge by 

trailing = trailing = <t> upper Slower 
edge edge 


From this starting value the circulation is convected along grid lines downstream of the 
trailing edge by 

r, + i>o 

The circulation is incorporated within the solution algorithm using the boundary 
condition 

A ^=0 


which is imposed numerically using 



upper 


~ (jk- i + r) 

Z k ~ Z k - 1 


and 



lower 


z k ~ z k - 1 
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To further demonstrate the accuracy of the ASP capability in producing solutions of 
similar accuracy to full potential solutions, two additional cases are presented for the 
NACA 0012 airfoil. The first case is a transonic application of the NACA 0012 airfoil at 

= 0.7 and a = 2°. The ASP pressure coefficient distribution is presented in Figure 10 
and two full potential solutions for this case reported by Wedan and South 25 are presented 
in Figure 11. Comparisons of the two figures indicate that the strength and location of 
the shock wave near 28% chord on the upper surface of the airfoil in the ASP result are in 
very good agreement with the FP results. Also, the ASP pressure distribution near the 
leading edge where stagnation occurs is smooth and well predicted in comparison with 
the FP pressure distributions. 

The second case is a high lift application of the NACA 0012 airfoil at M v = 0.3 and a 
= 10°. This case is generally regarded as being outside the range of applicability of small 
perturbation theory, and it is impossible computationally to obtain a converged solution 
with most small perturbation codes. The ASP pressure coefficient distribution is 
presented in Figure 12 and two full potential solutions reported by Wedan and South 25 are 
presented in Figure 13. Near the leading edge there is a rapid flow expansion about the 
upper surface and a small amount of reverse flow (u < 0) on the lower surface due to the 
relatively high angle of attack. The rapid expansion of the flow produces pressure 
coefficients close to the value of the sonic pressure coefficient near C P ~ -7 . The ASP 
pressure distribution is in very good agreement with the FP pressure distributions overall, 
including the prediction of the relatively large values of the pressure coefficients on the 
upper surface near the leading edge. 

Second-Order Accurate Supersonic Differencing 

Most small perturbation codes use first-order accurate differencing at supersonic 
points because of difficulty in constructing backward-biased second-order accurate 
difference stencils. The CAP-TSD code has an option for second-order accurate 
supersonic differencing but it does not work correctly because of a conceptual error in 
implementation. Therefore the first-order accurate differencing is used typically and 
shock waves are generally captured that are slightly weaker and located forward of their 
correct strength and location. In the ASP capability, a second-order accurate backward 
(upwind) calculation of the mass flux is achieved at cell interfaces for supersonic points. 
Similar to reconstruction of flow variables in advanced Euler and Navier-Stokes codes 
involving Riemann solvers, a flux limiter is required at the shock wave to prevent 
overshoots or oscillations in the solution because of the change in the solution there. 
Across the shock wave, although the potential (j) is continuous, the velocity u is 
discontinuous. The flux limiter allows the solution to be second-order accurate in smooth 
regions, and near the shock wave, the solution is formally first-order accurate. 

The effects of order of accuracy of supersonic differencing are demonstrated in 
Figure 14 for a case involving a strong shock wave. The case considered is a transonic 
application of the NACA 0012 airfoil at M v = 0.75 and a = 2°. Two sets of pressure 
coefficient distributions are presented corresponding to ASP calculations performed with 
first- and second-order accurate supersonic differencing. The two sets of pressures are in 
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very close agreement with each other with the only differences in the region of the shock 
wave on the upper surface of the airfoil. The strong shock in the first-order accurate 
solution is located near approximately 57% chord and the shock in the second-order 
solution is of similar strength and located near 60% chord. Figure 15 shows a full 
potential pressure distribution computed by Hafez 26 for this case using an accurate flux- 
biased method on a reasonably fine mesh. The FP distribution indicates that the strong 
shock wave is located correctly at approximately 60% chord on the upper surface of the 
airfoil. Thus, the ASP solution with second-order accurate supersonic differencing 
accurately predicts the shock location as well as the shock strength. 

Entropy and Vorticity Effects 

Shock-generated entropy effects are incorporated within small perturbation codes 
such as CAP-TSD by first using the Prandtl relation 17 (shock jump condition) to 
determine the velocity downstream of the shock wave from the upstream and sonic 
velocities. The upstream and downstream velocities are then used in the Rankine- 
Hugoniot (R-H) relation 17 to determine the change in entropy across the shock wave. 27 
The resulting change in entropy is subsequently used in a Clebsch formulation 28 to 
determine the vorticity downstream of the shock wave. 29,30 The vorticity modifies the 
calculation of the velocity field in the downstream region. This is a common procedure 
that has been used in various computer codes especially at the full-potential equation 
level. However, the approach as applied to the small perturbation equation, with any of 
the fluxes defined above including that of Williams, does not conserve mass. 

In contrast, the approach developed for the ASP capability conserves mass by using a 
ratio of the streamwise fluxes evaluated using the upstream and downstream velocities. 
Specifically, the downstream perturbation velocity is first computed using the Prandtl 
relation (shock jump relation) given by 


> [i +(<#■,) T 

\ _ |_ ' some J i 

where the subscripts 1 and 2 represent stations that are upstream and downstream of the 
shock wave, respectively. The change in entropy across the shock is then computed 
using 

As-cr-i){i-/ 1 [(^) l ]// l [(0 I ) 2 | 

For steady flows the entropy is held constant along gridlines downstream of shock waves. 
For unsteady flows the entropy is convected downstream using 

d . d 

— As + — As = 0 

dt dx 

The fluxes downstream of the shock are subsequently modified according to 
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which conserves mass across the shock wave by design. 

Similar to the CAP-TSD code, the ASP capability uses a Clebsch formulation 28 to 
compute the shock-generated vorticity. In brief, the streamwise velocities downstream of 
shocks are computed using 
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In the CAP-TSD code, when entropy and vorticity effects are included in the 
solution procedure no changes are made to the wake modeling because the first order 
term due to vorticity exactly cancels the first-order entropy term. 15 However, the second- 
order terms were found to be not insignificant, and thus, they were incorporated into the 
ASP solution procedure. Hence when entropy and vorticity are included, the circulation 
is convected using 
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where the subscripts “upper” and “lower” correspond to the values on the upper and 
lower surfaces of the trailing wake. 

The effects of entropy and vorticity modeling are demonstrated in Figure 16. The 
case considered is the NACA 0012 airfoil at M v = 0.75 and a = 2°, the same case 
presented in Figures 14 and 15. Here, three sets of ASP calculations were performed 
with the second-order accurate supersonic differencing including: (1) the isentropic 
capability (the same result shown in Figure 14 with the shock wave located at 60% 
chord), (2) the Rankine-Hugoniot shock jump relation to calculate the entropy change 
across the shock, and (3) the mass conserving calculation of the entropy change. When 
the R-H relation is used to calculate the entropy change the shock wave is weakened and 
located forward of the isentropic solution. However, when compared to the Euler 
calculation performed by Hafez 26 for this case shown in Figure 17, the R-H shock wave is 
not located far enough forward on the upper surface of the airfoil. Alternatively when the 
mass conserving calculation was used to determine the entropy change, a weaker shock 
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located near 46% chord was predicted, which is in very good agreement with the Euler 
result. Hence, the mass conserving calculation of the entropy change across the shock 
wave and the inclusion of the second-order terms in the trailing wake are essential 
elements in producing Euler-like accuracy with the ASP formulation. 

To further demonstrate the accuracy of the ASP capability in producing solutions of 
similar accuracy to Euler solutions, two additional cases are presented for the NACA 
0012 airfoil. The first case is a transonic application of the NACA 0012 airfoil at = 
0.6 and a = 5°. This is a difficult case for a small perturbation method because of the 
higher angle of attack. The ASP pressure coefficient distribution is presented in Figure 
18 and several solutions for this case reported by Klopfer and Nixon 31 are presented in 
Figure 19. The latter solutions include calculations involving FP, modified-FP, and Euler 
codes, the last of which is most relevant. Comparisons of the two figures indicate that the 
strength and location of the shock wave near 16% chord on the upper surface of the 
airfoil in the ASP result are in very good agreement with the Euler result. Also, the ASP 
pressure levels forward of the shock wave are well predicted in comparison with the 
Euler calculation. 

The second case is another transonic application of the NACA 0012 airfoil, this time 
at M v = 0.8 and a = 1.25°. The ASP pressure coefficient distribution is presented in 
Figure 20 and an Euler solution computed by Anderson, et al. 32 with the CFL3D code is 
presented in Figure 21. The case involves a strong shock wave on the upper surface of 
the airfoil at approximately 64% chord and a weak shock wave on the lower surface near 
34% chord. The difficulty of this case involves the difference in strength of the two 
shock waves. Comparisons of the two pressure distributions indicate that the ASP 
formulation including entropy and vorticity accurately predicts the strength and location 
of both of the shock waves, even though the shocks have disparate strengths. 

- Viscous Effects 

The dissipation integral method developed by Whitfield and coworkers 18 ’ 19 was 
implemented within the ASP3D code to model attached or shock-induced separated 
boundary layers. The dissipation integral approach represents a significant improvement 
over the Green’s lag entrainment method such as that in CAP-TSD 21 because the closure 
relations are derived by a detailed fitting of velocity profiles obtained through 
experiments for both attached and separated flows. Hence, the closure relations are 
continuous functions applicable to either attached or separated flows, and there is no need 
for any nonphysical interpolation between disparate sets of closure relations as done by 
Edwards. 21 Drela 20 has refined the dissipation integral method, wherein the closure 
relations depend not only upon the local Mach number and shape factor, but also on the 
momentum thickness Reynolds number. 

The ASP viscous method involves solving the unsteady boundary layer equations 
simultaneously with the outer potential flow solution so that no interaction law coupling 
the inner and outer solutions is required. The combined solution procedure involves an 
implicit block tridiagonal inversion for all of the cells along the surface and trailing wake. 
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Unlike the CAP-TSD viscous capability, 21 the ASP capability uses exact formulas for 
edge quantities and exact boundary conditions are used along surfaces and wakes. 
Smoothing of edge quantities and limiters are not required for stability, and no arbitrary 
or free parameters are necessary to tune the procedure. 

Specifically in the ASP capability, the exact edge quantities are used for the 
calculation of velocity u e , density p e , temperature T e , local Mach number M e , and the 
coefficient of viscosity u e given by 
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where S is Sutherland’s coefficient and T v is the freestream temperature. These equations 
are in contrast with the edge quantities contained within CAP-TSD wherein the density, 
temperature, and local Mach number relations were linearized as 
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For the ASP capability, the time-dependent integral boundary layer (IBL) and lag 
equations may be written as a system of equations in the form 
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where the independent variables are the momentum thickness 6, the incompressible shape 
parameter H k , and the square root of the shear stress coefficient c x . The matrix [ B ] has 
three eigenvalues that indicate the nature of the boundary layer. For attached boundary 
layers ( H k < H 0 ) all three eigenvalues are positive, as shown in Figure 22, indicating that 
all of the characteristics are in the downstream or positive direction. The shape 
parameter H 0 is defined as 20 

H 0 = 3 + 400 /Re 0 if Re e > 400 


H 0 = 4 if Re e < 400 

At separation ( H k = H 0 ) one of the eigenvalues becomes zero. And for separated flows 
( H k > H 0 ), one of the eigenvalues becomes negative, mathematically reflecting the fact 
that within the boundary layer there is some reverse flow (in the negative streamwise 
direction). 


where Re e - ^—^—Re 


For solution, the integral boundary layer and lag equations are first premultiplied by 
[A]' 1 to diagonalize the time term as 
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The resulting IBL system is then linearized in a simple way and cast into the so-called 
delta-form for implicit solution as 
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where the viscous residual is defined as 
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The system of equations is directly coupled with the ASP equation for the outer 
potential flow for simultaneous implicit solution using a block tridiagonal matrix 
inversion procedure. This is done for the cells adjacent to the upper and lower sides of 
the lifting surfaces and their wakes. All other cells do not involve the boundary layer, 
and hence, only the ASP equation is solved in those cells. Dissipation integral relations 
are used for turbulent closure, similar to those reported by Drela. 20 
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Surface Boundary Condition Including Entropy, Vorticity, and Viscous 
Effects 


In CAP-TSD the lowest order linear surface tangency condition is used including 
entropy, vorticity, and viscous effects. In ASP3D the complete mass-consistent surface 
boundary condition including entropy, vorticity, and viscous effects is given by 
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The surface is represented by 

b(x,t ) ± 6 *(x,t) - xa = 0 

where the ± is for the upper and lower surfaces, respectively, since the boundary layer 
displacement thickness <5* is defined as a positive quantity on both surfaces. 

Wake Boundary Condition Including Entropy, Vorticity, and Viscous Effects 

Viscous effects are included in the wake modeling in the ASP capability through the 
boundary condition 

A<fi z -±±(f,d'l+d; 
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across the wake where /i and g were defined previously. The boundary condition is 
incorporated numerically using 
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Three cases were considered to demonstrate the ASP viscous capability 
corresponding to an attached boundary layer, a shock-induced flow separation, and a 
naturally unsteady flow where the boundary layer periodically separates and reattaches. 
The attached flow case 33 is the RAE 2822 airfoil at = 0.676, a = -2.25°, and Re = 5.7 
x 10 6 , where the calculations were performed using a CFL number of thirty for five 
hundred time steps, and the results are presented in Figure 23. Figure 23(a) shows the 
convergence histories of the residuals for the outer potential flow R ^ and the inner IBF 
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equations R d , R H and R c on the upper surface of the airfoil. The residuals indicate that 

the solutions for both the outer ASP potential flow and the inner IBL indeed converge. 
The corresponding ASP pressure coefficient distribution is compared with the 
experimental pressure data 33 in Figure 23(b). These comparisons show very good 
agreement between the ASP pressures and the data, thus demonstrating the accuracy of 
the ASP viscous capability. 

The second case 34 corresponds to separated flow on the upper surface of the NACA 
0012 airfoil at M „ = 0.775, a = 2.05°, and Re = 10 7 . The calculations were performed 
with a CFL number of fifteen for a total of one thousand time steps and the results are 
presented in Figure 24. The CFL number for the separated flow case was half of that 
used in the attached flow case, although it is still a relatively large value. Figure 24(a) 
shows the convergence history in the more conventional form of the residual of the ASP 
potential equation. The convergence history indicates that the solution has converged 
approximately five orders of magnitude for the separated flow case. The corresponding 
ASP pressure coefficient distribution is shown in Figure 24(b), which is in good 
agreement with the experimental data, 34 especially in the strength and location of the 
upper surface shock wave. The shock induced separation along the upper surface of the 
airfoil ranges from the foot of the shock near 53% chord to the reattachment point at 
approximately 67% chord. 

The third case 33 that was considered to demonstrate the ASP viscous capability was 
the periodic flow separation about the RAE 2822 airfoil at M v = 0.75, a = 2.81°, and Re 
= 6.5 x 10 6 . The calculations were performed using a CFL number of twenty and the 
resulting lift and moment coefficient time histories are shown in Figure 25(a) for a total 
of three thousand time steps. The calculations were run much longer than the previous 
two cases to demonstrate clearly the periodicity of the coefficients as a function of time. 
In this case the flow on the upper surface of the airfoil periodically separates and later 
reattaches through each of the cycles shown. The instantaneous ASP pressure coefficient 
distribution at t = 3000 At is shown in Figure 24(b). The pressure distribution is smooth 
near the shock and the trailing edge, which are the two regions of strongest interaction 
between the boundary layer and the outer potential flow. Although not obvious from the 
figure, the flow is separated from the shock wave near 64% chord to the trailing edge. 


Introduction of the ASP3D Computer Program 

The ASP theory has been used as the basis for a new computer code called ASP3D 
(Advanced Small Perturbation - 3D), which involves either AF1- or AF2-type 
approximate factorization algorithms and a FAS (full approximation scheme) multigrid 
procedure for the solution of the ASP potential flow equation and associated boundary 
conditions. The ASP3D code can treat aircraft configurations involving multiple lifting 
surfaces with leading and trailing edge control surfaces and a fuselage. The new code is 
briefly described here and some representative results are presented in the following 
section. 
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Governing Equation 

The three-dimensional general small perturbation equation may be written in 
Cartesian coordinates as 

dt dx dy dz 

where the fluxes are defined as 
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The ASP3D code though solves the governing equation written in computational 
coordinates as 
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The various geometric quantities used in these equations such as the cell volume, areas of 
the cell interfaces, and direction cosines are defined elsewhere. 22 All of these quantities 
are computed using formulas that are exact. 


Entropy, Vorticity, and Viscous Effects 

Shock-generated entropy and vorticity effects are incorporated within the ASP3D 
code using the mass consistent modeling as described earlier. Viscous effects are 
included through a simultaneous implicit solution of the integral boundary layer and lag 
equations with the outer ASP potential flow as also described earlier. The ASP3D 
capability does not require the use of an interaction law because the equations are solved 
simultaneously. This is a stable approach that also does not require any smoothing or 
limiters. Turbulent closure is through the use of the dissipation integral relations of 
Drela, 20 which is applicable to attached and separated flows. 


Surface Boundary Condition 


In CAP-TSD the lowest order linear surface tangency condition is used including 
entropy, vorticity, and viscous effects. In ASP3D the complete mass-consistent surface 
boundary condition including entropy, vorticity, and viscous effects is given by 
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The surface is represented by 

b(x,y,t ) ± 6 *(x,t) - xa = 0 

where the ± is for the upper and lower surfaces, respectively, since the boundary layer 
displacement thickness 6* is defined as a positive quantity on both surfaces. 

Note that the spanwise slopes b y are included in the complete surface boundary 
condition for three-dimensional applications, which may be important for swept wings 
near the leading edge and also in the wing tip region. Furthermore, the streamwise and 
spanwise slopes are computed internal to ASP3D from the ordinates of the wing 
geometry by simple second-order accurate central finite differences, effectively the same 
as is done in advanced codes with body fitted meshes. This eliminates the differentiation 
of the spline of the wing geometry normally performed as a preprocessing step to 
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determine the airfoil slopes. It thus makes the input to the ASP3D code simpler than 
most other small perturbation codes such as CAP-TSD. 

Wake Boundary Condition 

Viscous effects are included in the wake modeling in the ASP capability through the 
boundary condition 

(y8*) x +«; 
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across the wake where / and g were defined previously. The boundary condition is 
incorporated numerically using 
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Shock Capturing 

The AF1 and AF2 algorithms of the ASP3D code use one of three approaches to 
model supersonic regions of the flow including: 1) Godunov, 35 2) Engquist-Osher 36 (E- 
O), and 3) Murman-Cole 37 (M-C) type-dependent switches. The first two switches satisfy 
an entropy inequality by design. This means that their use precludes the possibility of 
computing nonphysical expansion shocks that can occur with the Murman-Cole switch, 
especially in two-dimensional applications. With Godunov 35 switching, shocks are 
captured sharply, with only one interior state. Therefore, the Godunov switch is 
preferable over the E-0 and M-C switches, and it is the default switch within the ASP3D 
code. The E-0 switch is available for comparisons with CAP-TSD, and the M-C switch 
was included for completeness and historical reasons. The details of the ASP3D 
implementation of the shock capturing options are reported by Batina. 22 

Solution Algorithm 

The ASP3D code involves a modern cell-centered finite-volume discretization, which 
allows exact planform treatment including leading edges, wing tips, and control surface 
edges. The advanced small perturbation equation is solved using either an AF1 
approximate factorization algorithm or an AF2-type algorithm. The AF1 algorithm 
within the ASP3D code is the cell-centered finite-volume version of the AF1 algorithm 13 
developed for CAP-TSD, which may be written as 

hhh A0 = -*(0) 
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where the left hand side is an approximate factorization of a linearization of the 
governing equation, and R((p) is the residual. This equation is solved numerically in a 
three-sweep procedure wherein each operator is applied individually as 

L g A0 = -R 
= A0 


L^A<p = A (j) 


and the sweeps are performed independently in streamwise, spanwise, and vertical 
directions. The solution is advanced with a constant time step, which is selected for 
numerical stability. Specifically, the size of the time step is selected by trial and error so 
that shock waves in the flow do not move more than one cell per step. If that condition is 
violated, the solution becomes unstable, and the calculations diverge. This is the so- 
called moving shock instability inherent with AFl-type schemes as applied to the small 
perturbation equation. 

The AF2 algorithm within the ASP3D code is described in general by 

hhh A0 = -/e(0) 


The general form of the AF2 algorithm appears similar to that of the AF1 algorithm but 
there are two fundamental differences. First, the operators are applied in the reverse 
order (vertical sweep and spanwise sweep performed first, with the streamwise sweep 
performed last) and a streamwise implicit temporal damping term is included on the right 
hand side of the equation which is indicated below. The AF2 algorithm is defined by 
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In the first sweep, the temporal damping term is treated implicitly since it involves 
information from the second intermediate sweep but at the previous grid plane i - 1, and 
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a = — + w 
At 


^Ka + ^b 


The parameter ip is set equal to zero for steady applications and set equal to unity for 
unsteady applications. The parameter co is an over-relaxation factor. The algorithm 
involves both physical (At) and computational (At) time steps. The computational step is 
selected for rapid convergence and it is usually assigned a large value. Hence the 
solution can be advanced using either a constant computational time step or a constant 
CFL number. The physical time step is used only for unsteady applications, and it is 
selected based on the physics of the problem being considered. 

The AF2 approach is the preferred algorithm. The AF2 scheme is more robust than 
the AF1 scheme since it cannot fail due to the moving shock instability by design. The 
temporal damping term eliminates the moving shock instability. Hence convergence to 
steady state is achieved with large time steps or large CFL numbers, independent of the 
moving shock instability that plagues the AF1 algorithm. The AF2 scheme may also be 
used for unsteady calculations when local iteration is applied. Therefore, the physical 
time step may be selected based on the physics of the problem being solved, rather than 
on numerical stability considerations. For either steady or unsteady applications, the 
temporal damping term vanishes in the steady-state limit or within the convergence 
process of performing subiterations for unsteady calculations. 

Multigrid Implementation 

The AF1 algorithm is not appropriate for multigrid implementation because the rapid 
rate of convergence of the multigrid procedure would result in the moving shock 
instability. The AF2 algorithm of the ASP3D code is amenable to multigrid 
implementation because the temporal damping term prevents the moving shock 
instability from occurring. Therefore the multigrid procedure 38 may be used for steady or 
iterative unsteady applications with ASP3D. The procedure is implemented in FAS (Full 
Approximation Scheme) form, 39 which is applicable to nonlinear governing equations. 
As with any multigrid procedure, high-frequency errors are damped on the fine mesh and 
low-frequency errors are damped on the coarser meshes. Either Y- or W-cycles may be 
selected and full multigrid is available. W-cycles have been found to be slightly more 
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effective for solution convergence, although W-cycles are slightly more expensive than 
V-cycles. 40 

Specifically, the calculation on the finest mesh (N) smoothes high frequency errors 
using the AF2 algorithm written in general form 

Subsequent calculations on the next coarser mesh (AM) involve fine mesh residual 
injection 

^v-i(^v-i) = ^N-i ~In {^n) + ^n-Jn {Qn) 


where the I functions are restriction operators for transferring the velocity potential and 
residual from the fine mesh to the coarse mesh. The procedure is generally repeated for 
additional coarse meshes ( N-2 , N-3, . . .). 


The restriction operator 40 for potential is a volume weighted operator defined by 




and the restriction operator 40 for residuals is defined by 
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where the summations above are taken over all of the fine mesh cells that make up the 
coarse mesh cell. The interpolation operator for the velocity potential that is required to 
transfer information back from coarse to fine mesh is a trilinear interpolation operator 
defined symbolically by 

&N &N + 


ASP3D Results and Discussion 

Representative ASP3D results are presented briefly to demonstrate application of the 
ASP theory to more practical problems. Two cases were considered including: (1) 
inviscid and viscous calculations for the F-5 fighter wing 41 and (2) viscous calculations 
for the NACA RM L51F07 wing-fuselage configuration. 42 All of the results are 
compared with experimental data to assess the accuracy of the ASP3D computations. 

- F-5 Fighter Wing Results 

The F-5 fighter wing 41 has a leading edge sweep of 31.92 degrees, an aspect ratio of 
1.578, and a taper ratio of 0.283. The wing has a NACA 65A004.8 airfoil section that 
has a drooped leading edge. The wing was modeled for ASP3D applications using a 97 x 
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25 x 65 finite volume mesh as shown in Figure 26. A near field view of the 97 x 25 
planform mesh is shown in Figure 26(a) and a view of the 97 x 65 root sectional mesh is 
shown in Figure 26(b). The meshes were generated using a program that was created to 
construct meshes specifically for use by the ASP3D code, and as such, it can generate 
meshes for aircraft configurations with multiple lifting surfaces, with multiple leading 
and trailing edge controls, and a fuselage. The mesh generator automatically creates the 
input file for ASP3D, which makes it user friendly. 

Inviscid and viscous calculations were performed both using the entropy and vorticity 
effects for the F-5 wing at M x = 0.946, a = -0.004°, and Re = 5.89 x 10 6 . These results 
are presented in Figure 27, with the inviscid pressure distributions shown in Figure 27(a) 
and the viscous pressure distributions shown in Figure 27(b). The ASP3D pressure 
distributions were interpolated to the eight span stations where the NLR experimental 
data 41 were measured ranging from r| = 0.174 in the inboard region of the wing to r| = 
0.939 near the wing tip. The inviscid pressure distributions (Figure 27(a)) show 
generally good agreement with the experimental pressure data except in the shock region 
on the upper and lower surfaces of the wing. This is expected because of the neglect of 
viscous effects. The viscous pressure distributions (Figure 27(b)) are similar to the 
inviscid pressures and thus also agree well with the data, although the shock waves are 
slightly weaker and located more forward, and hence, they are in better agreement with 
the experimental shock results. 

NACA RM L51F07 Wing-Fuselage Results 

The geometry of the NACA RM L51F07 wing-fuselage configuration 42 is defined in 
Figure 28. The wing has a leading edge sweep of 46.76 degrees (quarter chord sweep of 
45 degrees), an aspect ratio of 4.0, and a taper ratio of 0.6. The wing has a NACA 
65A006 airfoil section. The fuselage is an axisymmetric body with fineness ratio of 12. 
The model was sting mounted for testing and pressures were measured along five chords 
of the wing as well as on the fuselage. The NACA RM L51F07 configuration was 
modeled for ASP3D applications using a 129 x 29 x 75 finite volume mesh. A near field 
view of the 129 x 29 planform mesh is shown in Figure 29. The mesh was generated in 
the streamwise direction by clustering points near the wing, and there is a slight 
clustering near the fuselage nose and where the fuselage is attached to the sting. The 
mesh was generated in the spanwise direction by clustering cells near the fuselage and 
near the wing tip. 

Viscous calculations were performed using the entropy and vorticity effects for the 
NACA RM L51F07 configuration at M x = 0.93, a = 2.0°, and Re = 10 7 . These results are 
presented in Figure 30 including comparisons with experiment data, 42 with the wing 
pressure comparisons shown in Figure 30(a) and the fuselage pressure comparisons 
shown in Figure 30(b). The ASP3D pressure distributions on the wing were interpolated 
to the five span stations where experimental data 42 were measured ranging from r| = 0.2 
near the wing fuselage junction to r| = 0.95 near the wing tip. The ASP3D pressure 
distributions on the fuselage were directly computed along the fuselage centerline for 
comparison with the experimental data measured at r| = 0.0. The wing pressure 
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comparisons of Figure 30(a) show good overall agreement in pressure levels and in the 
prediction of the strength and location of the upper surface shock wave. The largest 
differences between computation and experiment are in the inboard region at r| = 0.2, 
which may be the case because there is no boundary layer on the fuselage (only on the 
wing). The fuselage pressure comparisons of Figure 30(b) also show good agreement 
although the calculated shock wave is a little sharper than the wing shock wave due to the 
neglect of viscous effects on the fuselage. 

Conclusions 

A new advanced small perturbation (ASP) potential flow theory was developed by 
first determining the essential elements required to produce solutions as accurate as a full 
potential code with the small perturbation approach on a Cartesian grid. This level of 
accuracy was obtained by using a higher-order streamwise mass flux and a mass 
conserving surface boundary condition. Subsequent applications with these 
improvements showed very good agreement with all of the full potential cases 
considered. This was true even for cases at freestream Mach numbers as low as 0.3 and 
angles of attack as high as ten degrees, conditions that are normally considered to be 
outside the range of applicability of classical small perturbation theories. 

The ASP theory was further developed to produce Euler-like solutions by 
incorporating a mass conserving calculation of the entropy jump across of the shock 
wave based on the higher-order streamwise mass flux of the ASP theory. Second-order 
terms in the trailing wake boundary condition were found to be not insignificant for 
lifting cases and were also incorporated. Subsequent calculations with these 
modifications showed very good agreement with all of the Euler cases that were 
considered. This was true for cases involving strong shock waves and multiple shock 
waves of disparate strengths. 

Viscous effects were incorporated within the ASP formulation by coupling an 
integral boundary layer procedure to the outer inviscid flow calculation. The dissipation 
integral method was implemented to model attached or shock-induced separated 
boundary layers. The capability involves solving the unsteady boundary layer and lag 
equations simultaneously with the outer potential flow solution, so that no interaction law 
coupling the inner and outer solutions is required. The combined solution procedure 
involves an implicit block tridiagonal inversion for all of the cells along the surface and 
its trailing wake. Exact formulas for edge quantities and exact boundary conditions are 
used along surfaces and wakes. Smoothing of edge quantities and limiters are not 
required for stability, and no arbitrary or free parameters are necessary to tune the 
procedure. Calculations performed with the complete ASP capability including entropy, 
vorticity, and viscous effects showed good agreement with experimental data. 

The ASP theory was shown to be mathematically more appropriate and 
computationally more accurate than the classical TSP theories. Consequently, the new 
theory has been used as the basis for a new computer code called ASP3D (Advanced 
Small Perturbation - 3D), which involves either AF1- or AF2-type approximate 
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factorization algorithms and a FAS (full approximation scheme) multigrid procedure for 
the solution of the ASP potential flow equation. The ASP3D code can treat aircraft 
configurations involving multiple lifting surfaces with leading and trailing edge control 
surfaces and a fuselage. The new code was briefly described and some representative 
results were presented. Detailed descriptions of the ASP3D capability and 
comprehensive results and comparisons are reported in a companion publication. 22 
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Table 1 - Mass flux, derivative of mass flux, and sonic velocity for various small perturbation theories. 
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Table 2 - Speed of sound, velocity, and stagnation quantities for various small perturbation theories. 



TSP theory with 
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U = l + <f>x 


(a) Stream wise mass flux/,. 



u = l + <f> x 

(b) Derivative of stream wise mass flux df { /d^ x . 

Figure 1 - Mass flux quantities as functions of stream wise velocity u for various 
small perturbation theories at M„ = 0.72. 
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u = i + 4> x 

(a) Stream wise mass flux/,. 



U = 1 + <t>X 

(b) Derivative of stream wise mass flux df,/df x . 

Figure 2 - ASP mass flux quantities as functions of stream wise velocity u for a range 

of values of freestream Mach number. 
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a/i/e* 



» XU 

1.30 1.35 1.40 1.45 1.50 


u = 1 + 0 * 

Figure 3 - Sonic point ( df,ldf x = 0) as a function of streamwise velocity u for various 

small perturbation theories at = 0.72. 


34 





Figure 4 - Cartesian finite volume mesh (257 x 129) with fifty chord extent for 
airfoil applications using the various small perturbation theories. 



(a) Leading edge region. 


(b) Trailing edge region. 


Figure 5 - Near field views of the finite volume mesh for airfoil applications using 
the various small perturbation theories. 
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Figure 6 - Effects of small perturbation streamwise flux on pressure coefficient 
distribution for the NACA 64A410 airfoil at M„ = 0.72 and a = 0°. 



shock 

location 


Figure 7 - Full potential pressure coefficient distribution computed by Jameson 24 
for the NACA 64A410 airfoil at M. = 0.72 and o = 0°. 
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Figure 8 - Effects of small perturbation surface boundary condition on pressure 
coefficient distribution for the NACA 64A410 airfoil at M„ = 0.72 and a = 0°. 
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Figure 9 - Full potential pressure coefficient distribution computed by Jameson 24 
for the NACA 64A410 airfoil at = 0.72 and o = 0° (same as Figure 7). 
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Figure 10 - ASP pressure coefficient distribution for the NACA 0012 airfoil at 

M m = 0.7 and a = 2°. 



Figure 11 - Full potential pressure coefficient distributions reported by Wedan and 
South 25 for the NACA 0012 airfoil at M m = 0.7 and a = 2°. 
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Figure 12 - ASP pressure coefficient distribution for the NACA 0012 airfoil at 

M„ = 0.3 and a = 10°. 



Figure 13 - Full potential pressure coefficient distributions reported by Wedan and 
South 25 for the NACA 0012 airfoil at = 0.3 and a = 10°. 
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Figure 14 - Effects of order of accuracy of supersonic differencing on pressure 
coefficient distribution for the NACA 0012 airfoil at M m = 0.75 and a = 2°. 



Figure 15 - Full potential pressure coefficient distribution computed by Hafez 26 for 
the NACA 0012 airfoil at = 0.75 and a = 2°. 
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Figure 16 - Effects of small perturbation entropy and vorticity modeling on pressure 
coefficient distribution for the NACA 0012 airfoil at M m = 0.75 and a = 2°. 

-a.ol 

shock 

-i.»- location 



Figure 17 - Euler pressure coefficient distribution computed by Hafez 26 for the 
NACA 0012 airfoil at M m = 0.75 and a = 2°. 


41 




Figure 18 - Pressure coefficient distribution for the NACA 0012 airfoil at M„ = 0.6 
and a = 5° computed using ASP theory including entropy and vorticity. 



Figure 19 - Euler and FP pressure coefficient distributions reported by Klopfer and 
Nixon 31 for the NACA 0012 airfoil at = 0.6 and a = 5°. 
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Figure 20 - Pressure coefficient distribution for the NACA 0012 airfoil at M„ = 0.8 
and a = 1.25° computed using ASP theory including entropy and vorticity. 



x/c 

Figure 21 - Euler (CFL3D) pressure coefficient distribution computed by Anderson, 
et al. 32 for the NACA 0012 airfoil at M m = 0.8 and a = 1.25°. 
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10 ° 10 1 10 2 

H k 


Figure 22 - Normalized eigenvalues of the integral boundary layer and lag 
equations as functions of the incompressible shape parameter H k for Re 0 = 400 

and = 1.0. 
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(a) Convergence histories for outer potential flow and inner IBL solutions. 



x/c 

(b) ASP pressure coefficient comparisons with experimental data. 33 

Figure 23 - ASP viscous calculations with CFL = 30 for the RAE 2822 airfoil at 
M x = 0.676, a = -2.25°, and Re = 5.7 x 10 6 (attached flow case). 
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iterations 


(a) Convergence history for the outer potential flow solution. 



x/c 

(b) ASP pressure coefficient comparisons with experimental data. 34 

Figure 24 - ASP viscous calculations with CFL = 15 for the NACA 0012 airfoil at 
M x = 0.775, a = 2.05°, and Re = 10 7 (separated flow case). 
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(a) Lift and moment coefficient time histories with CFL = 20. 

- 1.6 
- 1.2 
- 0.8 
- 0.4 
C p 0.0 
0.4 
0.8 
1.2 

- 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 

x/c 

(b) Instantaneous pressure coefficient distribution at t = 3000At. 

Figure 25 - ASP viscous calculations for the RAE 2822 airfoil at M„ = 0.75, 
a = 2.81°, and Re = 6.5 x 10 6 (periodic flow case). 
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(b) Near field view of 97 x 65 root sectional mesh. 
Figure 26 - Finite volume meshes for the F-5 fighter wing. 41 
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ASP3D inviscid solution 

o NLR data - upper surface 
n NLR data - lower surface 





(b) Viscous calculation. 


Figure 27 - ASP3D pressure coefficient comparisons with NLR experimental data 41 
for the F-5 fighter wing at M x = 0.946, a = -0.004°, and Re = 5.89 x 10 6 . 
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Wing definition: 

- sweep =46.76 degrees 

- aspect ratio = 4.0 

- taper ratio = 0.6 

- NACA 65 A006 airfoil 


Fuselage definition: 

- axisymmetric body 

- fineness ratio = 12 

- sting mounted 


T] - 0 .0 



stations where 
pressure data 
were 

measured 


Figure 28 - Geometrical definition of the NACA RM L51F07 wing-fuselage 

configuration. 42 



Figure 29 - Near field view of 129 x 29 planform mesh of the 129 x 29 x 75 total 
finite volume mesh for the NACA RM L51F07 wing-fuselage configuration. 
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(a) Comparisons on the wing. 



NACA RM L51F07 
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(b) Comparisons on the fuselage. 


Figure 30 - ASP3D pressure coefficient comparisons with experimental data for the 
NACA RM L51F07 wing-fuselage configuration at = 0.93, a = 2.0°, and Re = 10 7 . 
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